﻿using System;

// http://www.cnblogs.com/GIS_zhou/articles/1525772.html
// http://wiki.openstreetmap.org/wiki/Mercator#C.23
// http://zh.wikipedia.org/wiki/%E9%BA%A5%E5%8D%A1%E6%89%98%E6%8A%95%E5%BD%B1%E6%B3%95

namespace DegreeDistance
{

    public class Pixel
    {
        public double X { get; set; }
        public double Y { get; set; }
    }

    public class GeoCoord
    {
        public double Lon { get; set; }
        public double Lat { get; set; }
    }

    public static class Mercator
    {
        private const double RMajor = 6378137.0;
        private const double RMinor = 6356752.3142;
        private const double Ratio = RMinor / RMajor;
        private static readonly double Eccent = Math.Sqrt(1.0 - (Ratio * Ratio));
        private static readonly double Com = 0.5 * Eccent;

        private const double Deg2Rad = Math.PI / 180.0;
        private const double Rad2Deg = 180.0 / Math.PI;
        private const double Pi2 = Math.PI / 2.0;

        public static Pixel ToPixel(double lon, double lat)
        {
            return new Pixel { X = LonToX(lon), Y = LatToY(lat) };
        }

        public static GeoCoord ToGeoCoord(double x, double y)
        {
            return new GeoCoord { Lon = XToLon(x), Lat = YToLat(y) };
        }

        public static double LonToX(double lon)
        {
            return RMajor * DegToRad(lon);
        }

        public static double LatToY(double lat)
        {
            lat = Math.Min(89.5, Math.Max(lat, -89.5));
            double phi = DegToRad(lat);
            double sinphi = Math.Sin(phi);
            double con = Eccent * sinphi;
            con = Math.Pow(((1.0 - con) / (1.0 + con)), Com);
            double ts = Math.Tan(0.5 * ((Math.PI * 0.5) - phi)) / con;
            return 0 - RMajor * Math.Log(ts);
        }

        public static double XToLon(double x)
        {
            return RadToDeg(x) / RMajor;
        }

        public static double YToLat(double y)
        {
            double ts = Math.Exp(-y / RMajor);
            double phi = Pi2 - 2 * Math.Atan(ts);
            double dphi = 1.0;
            int i = 0;
            while ((Math.Abs(dphi) > 0.000000001) && (i < 15))
            {
                double con = Eccent * Math.Sin(phi);
                dphi = Pi2 - 2 * Math.Atan(ts * Math.Pow((1.0 - con) / (1.0 + con), Com)) - phi;
                phi += dphi;
                i++;
            }
            return RadToDeg(phi);
        }

        private static double RadToDeg(double rad)
        {
            return rad * Rad2Deg;
        }

        private static double DegToRad(double deg)
        {
            return deg * Deg2Rad;
        }
    }

    public static class Util
    {
        /// <summary>
        /// 转换为精度5米的整形数
        /// </summary>
        /// <param name="a"></param>
        /// <returns></returns>
        public static int ConvertPrecision5(double a)
        {
            return Convert.ToInt32(Math.Round(a / 5.0d, 0));
        }
    }

}
